#pragma once
#include "DXUT.h"
typedef union{
	double dp;
	unsigned short sh[4];
} Double;

typedef union{
	double dp;
	unsigned short sh[2];
} Float;


#define MSW 0 /* 0 if the left-most 16-bit short is sh[0] */
/* 3 if the left-most 16-bit short is sh[3] */

static unsigned short mask[16] = { /* bit masks for bits 0 - 15 */
	0x0001, 0x0002, 0x0004, 0x0008, 0x0010, 0x0020, 0x0040, 0x0080,
	0x0100, 0x0200, 0x0400, 0x0800, 0x1000, 0x2000, 0x4000, 0x8000};

double ulpsoft(double x)
{
	double ulp;
	int exp;
	frexp(x,&exp);
	ulp=ldexp(0.5,exp-52);
	return ulp;
}

float ulpfsoft(float x)
{
	float ulp;
	int exp;
	frexp(x,&exp);
	ulp=ldexp(0.5f,exp-23);
	return ulp;
}

	
double ulp(double x)
{
	Double U,							/* ulp of x */
		X;								/* working copy of x */
	int bit,							/* position of bit e-1 in 16-bit word */
		e1,								/* biased exponent - 1 */
		word;							/* index of 16-bit word containing bit e-1 */
	X.dp = x;
	X.sh[MSW] &= 0x7ff0;				/* isolate exponent in 16-bit word */
										/* X.sh[0] now holds the exponent in bits 14-4 */
	U.dp = 0.0;							/* initialize exponent and mantissa to 0 */
	if (X.sh[MSW] > 0x0340)				/* ulp is normalized number */
		U.sh[MSW] = X.sh[MSW]-0x0340;	/* set exponent to e-52 */
										/* the value 0x0340 is 52 left-shifted 4 bits,	i.e. 0x0340 = 832 = 52<<4 */

	else {								/* ulp is denormalized number */
		e1 = (X.sh[MSW]>>4) - 1;		/* biased exponent - 1 */
		word = e1>>4;					/* find 16-bit word containing bit e-1 */
		if (MSW == 0) word = 3 - word;	/* compensate for word ordering */
		bit = e1%16;					/* find the bit position in this word */
		U.sh[word] |= mask[bit];		/* set the bit to 1 */
	}
	return U.dp;						/* return ulp */
}


//float ulpf(float x)
//{
//	Float U,X;
//	int bit, e1, word;
//	X.dp = x;
//	X.sh[MSW] &= ;
//}
